home *** CD-ROM | disk | FTP | other *** search
/ AGA Toolkit '97 / The AGA Toolkit '97.iso / miscellaneous / science / maths / calc / source / qmod.c < prev    next >
Encoding:
C/C++ Source or Header  |  1996-09-07  |  10.7 KB  |  484 lines
  1. /*
  2.  * Copyright (c) 1993 David I. Bell
  3.  * Permission is granted to use, distribute, or modify this source,
  4.  * provided that this copyright notice remains intact.
  5.  *
  6.  * Modular arithmetic routines for normal numbers, and also using
  7.  * the faster REDC algorithm.
  8.  */
  9.  
  10. #include "qmath.h"
  11.  
  12.  
  13. /*
  14.  * Structure used for caching REDC information.
  15.  */
  16. typedef struct    {
  17.     NUMBER    *num;        /* modulus being cached */
  18.     REDC    *redc;        /* REDC information for modulus */
  19.     long    age;        /* age counter for reallocation */
  20. } REDC_CACHE;
  21.  
  22.  
  23. static long redc_age;            /* current age counter */
  24. static REDC_CACHE redc_cache[MAXREDC];    /* cached REDC info */
  25.  
  26.  
  27. static REDC *qfindredc MATH_PROTO((NUMBER *q));
  28.  
  29.  
  30. /*
  31.  * Return the remainder when one number is divided by another.
  32.  * The second argument cannot be negative.  The result is normalized
  33.  * to lie in the range 0 to q2, and works for fractional values as:
  34.  *    qmod(q1, q2) = q1 - (int(q1 / q2) * q2).
  35.  */
  36. NUMBER *
  37. qmod(q1, q2)
  38.     register NUMBER *q1, *q2;
  39. {
  40.     ZVALUE res;
  41.     NUMBER *q, *tmp;
  42.  
  43.     if (qisneg(q2) || qiszero(q2))
  44.         math_error("Non-positive modulus");
  45.     if (qisint(q1) && qisint(q2)) {        /* easy case */
  46.         zmod(q1->num, q2->num, &res);
  47.         if (ziszero(res)) {
  48.             zfree(res);
  49.             return qlink(&_qzero_);
  50.         }
  51.         if (zisone(res)) {
  52.             zfree(res);
  53.             return qlink(&_qone_);
  54.         }
  55.         q = qalloc();
  56.         q->num = res;
  57.         return q;
  58.     }
  59.     q = qquo(q1, q2);
  60.     tmp = qmul(q, q2);
  61.     qfree(q);
  62.     q = qsub(q1, tmp);
  63.     qfree(tmp);
  64.     if (qisneg(q)) {
  65.         tmp = qadd(q2, q);
  66.         qfree(q);
  67.         q = tmp;
  68.     }
  69.     return q;
  70. }
  71.  
  72.  
  73. /*
  74.  * Given two numbers (a and b), calculate the quotient (c) and remainder (d)
  75.  * when a is divided by b.  This is defined so 0 <= d < b, and c is integral,
  76.  * and a = b * c + d.  The results are returned indirectly through pointers.
  77.  * This works for integers or fractions.  Returns whether or not there is a
  78.  * remainder.  Examples:
  79.  *    qquomod(11, 4, &x, &y) sets x to 2, y to 3, and returns TRUE.
  80.  *    qquomod(-7, 3, &x, &y) sets x to -3, y to 2, and returns TRUE.
  81.  */
  82. BOOL
  83. qquomod(q1, q2, retqdiv, retqmod)
  84.     NUMBER *q1, *q2;        /* numbers to do quotient with */
  85.     NUMBER **retqdiv;        /* returned quotient */
  86.     NUMBER **retqmod;        /* returned modulo */
  87. {
  88.     NUMBER *qq, *qm, *tmp;
  89.  
  90.     if (qisneg(q2) || qiszero(q2))
  91.         math_error("Non-positive modulus");
  92.  
  93.     if (qisint(q1) && qisint(q2)) {        /* integer case */
  94.         qq = qalloc();
  95.         qm = qalloc();
  96.         zdiv(q1->num, q2->num, &qq->num, &qm->num);
  97.         if (!qisneg(q1) || qiszero(qm)) {
  98.             *retqdiv = qq;
  99.             *retqmod = qm;
  100.             return !qiszero(qm);
  101.         }
  102.  
  103.         /*
  104.          * Need to fix up negative results.
  105.          */
  106.         tmp = qdec(qq);
  107.         qfree(qq);
  108.         qq = tmp;
  109.         tmp = qsub(q2, qm);
  110.         qfree(qm);
  111.         qm = tmp;
  112.         *retqdiv = qq;
  113.         *retqmod = qm;
  114.         return TRUE;
  115.     }
  116.  
  117.     /*
  118.      * Fractional case.
  119.      */
  120.     qq = qquo(q1, q2);
  121.     tmp = qmul(qq, q2);
  122.     qm = qsub(q1, tmp);
  123.     qfree(tmp);
  124.     if (qisneg(qm)) {
  125.         tmp = qadd(qm, q2);
  126.         qfree(qm);
  127.         qm = tmp;
  128.         tmp = qdec(qq);
  129.         qfree(qq);
  130.         qq = tmp;
  131.     }
  132.     *retqdiv = qq;
  133.     *retqmod = qm;
  134.     return !qiszero(qm);
  135. }
  136.  
  137.  
  138. #if 0
  139. /*
  140.  * Return the product of two integers modulo a third integer.
  141.  * The result is in the range 0 to q3 - 1 inclusive.
  142.  *    q4 = (q1 * q2) mod q3.
  143.  */
  144. NUMBER *
  145. qmulmod(q1, q2, q3)
  146.     NUMBER *q1, *q2, *q3;
  147. {
  148.     NUMBER *q;
  149.  
  150.     if (qisneg(q3) || qiszero(q3))
  151.         math_error("Non-positive modulus");
  152.     if (qisfrac(q1) || qisfrac(q2) || qisfrac(q3))
  153.         math_error("Non-integers for qmulmod");
  154.     if (qiszero(q1) || qiszero(q2) || qisunit(q3))
  155.         return qlink(&_qzero_);
  156.     q = qalloc();
  157.     zmulmod(q1->num, q2->num, q3->num, &q->num);
  158.     return q;
  159. }
  160.  
  161.  
  162. /*
  163.  * Return the square of an integer modulo another integer.
  164.  * The result is in the range 0 to q2 - 1 inclusive.
  165.  *    q2 = (q1^2) mod q2.
  166.  */
  167. NUMBER *
  168. qsquaremod(q1, q2)
  169.     NUMBER *q1, *q2;
  170. {
  171.     NUMBER *q;
  172.  
  173.     if (qisneg(q2) || qiszero(q2))
  174.         math_error("Non-positive modulus");
  175.     if (qisfrac(q1) || qisfrac(q2))
  176.         math_error("Non-integers for qsquaremod");
  177.     if (qiszero(q1) || qisunit(q2))
  178.         return qlink(&_qzero_);
  179.     if (qisunit(q1))
  180.         return qlink(&_qone_);
  181.     q = qalloc();
  182.     zsquaremod(q1->num, q2->num, &q->num);
  183.     return q;
  184. }
  185.  
  186.  
  187. /*
  188.  * Return the sum of two integers modulo a third integer.
  189.  * The result is in the range 0 to q3 - 1 inclusive.
  190.  *    q4 = (q1 + q2) mod q3.
  191.  */
  192. NUMBER *
  193. qaddmod(q1, q2, q3)
  194.     NUMBER *q1, *q2, *q3;
  195. {
  196.     NUMBER *q;
  197.  
  198.     if (qisneg(q3) || qiszero(q3))
  199.         math_error("Non-positive modulus");
  200.     if (qisfrac(q1) || qisfrac(q2) || qisfrac(q3))
  201.         math_error("Non-integers for qaddmod");
  202.     q = qalloc();
  203.     zaddmod(q1->num, q2->num, q3->num, &q->num);
  204.     return q;
  205. }
  206.  
  207.  
  208. /*
  209.  * Return the difference of two integers modulo a third integer.
  210.  * The result is in the range 0 to q3 - 1 inclusive.
  211.  *    q4 = (q1 - q2) mod q3.
  212.  */
  213. NUMBER *
  214. qsubmod(q1, q2, q3)
  215.     NUMBER *q1, *q2, *q3;
  216. {
  217.     NUMBER *q;
  218.  
  219.     if (qisneg(q3) || qiszero(q3))
  220.         math_error("Non-positive modulus");
  221.     if (qisfrac(q1) || qisfrac(q2) || qisfrac(q3))
  222.         math_error("Non-integers for qsubmod");
  223.     if (q1 == q2)
  224.         return qlink(&_qzero_);
  225.     q = qalloc();
  226.     zsubmod(q1->num, q2->num, q3->num, &q->num);
  227.     return q;
  228. }
  229.  
  230.  
  231. /*
  232.  * Return the negative of an integer modulo another integer.
  233.  * The result is in the range 0 to q2 - 1 inclusive.
  234.  *    q2 = (-q1) mod q2.
  235.  */
  236. NUMBER *
  237. qnegmod(q1, q2)
  238.     NUMBER *q1, *q2;
  239. {
  240.     NUMBER *q;
  241.  
  242.     if (qisneg(q2) || qiszero(q2))
  243.         math_error("Non-positive modulus");
  244.     if (qisfrac(q1) || qisfrac(q2))
  245.         math_error("Non-integers for qnegmod");
  246.     if (qiszero(q1) || qisunit(q2))
  247.         return qlink(&_qzero_);
  248.     q = qalloc();
  249.     znegmod(q1->num, q2->num, &q->num);
  250.     return q;
  251. }
  252. #endif
  253.  
  254.  
  255. /*
  256.  * Return the integer congruent to an integer whose absolute value is smallest.
  257.  * This is a unique integer in the range int((q2-1)/2 to int(q2/2), inclusive.
  258.  * For example, for a modulus of 7, returned values are [-3, 3], and for a
  259.  * modulus of 8, returned values are [-3, 4].
  260.  */
  261. NUMBER *
  262. qminmod(q1, q2)
  263.     NUMBER *q1, *q2;
  264. {
  265.     NUMBER *q;
  266.  
  267.     if (qisneg(q2) || qiszero(q2))
  268.         math_error("Non-positive modulus");
  269.     if (qisfrac(q1) || qisfrac(q2))
  270.         math_error("Non-integers for qminmod");
  271.     if (qiszero(q1) || (q1->num.len < q2->num.len - 1))
  272.         return qlink(q1);
  273.     q = qalloc();
  274.     zminmod(q1->num, q2->num, &q->num);
  275.     return q;
  276. }
  277.  
  278.  
  279. /*
  280.  * Return whether or not two integers are congruent modulo a third integer.
  281.  * Returns TRUE if the numbers are not congruent, and FALSE if they are.
  282.  */
  283. BOOL
  284. qcmpmod(q1, q2, q3)
  285.     NUMBER *q1, *q2, *q3;
  286. {
  287.     if (qisneg(q3) || qiszero(q3))
  288.         math_error("Non-positive modulus");
  289.     if (qisfrac(q1) || qisfrac(q2) || qisfrac(q3))
  290.         math_error("Non-integers for qcmpmod");
  291.     if (q1 == q2)
  292.         return FALSE;
  293.     return zcmpmod(q1->num, q2->num, q3->num);
  294. }
  295.  
  296.  
  297. /*
  298.  * Convert an integer into REDC format for use in faster modular arithmetic.
  299.  * The number can be negative or out of modulus range.
  300.  */
  301. NUMBER *
  302. qredcin(q1, q2)
  303.     NUMBER *q1;        /* number to convert into REDC format */
  304.     NUMBER *q2;        /* modulus */
  305. {
  306.     REDC *rp;        /* REDC information */
  307.     NUMBER *r;        /* result */
  308.  
  309.     if (qisfrac(q1))
  310.         math_error("Non-integer for qredcin");
  311.     rp = qfindredc(q2);
  312.     if (qiszero(q1))
  313.         return qlink(&_qzero_);
  314.     r = qalloc();
  315.     zredcencode(rp, q1->num, &r->num);
  316.     return r;
  317. }
  318.  
  319.  
  320. /*
  321.  * Convert a REDC format number back into a normal integer.
  322.  * The resulting number is in the range 0 to the modulus - 1.
  323.  */
  324. NUMBER *
  325. qredcout(q1, q2)
  326.     NUMBER *q1;        /* number to convert out of REDC format */
  327.     NUMBER *q2;        /* modulus */
  328. {
  329.     REDC *rp;        /* REDC information */
  330.     NUMBER *r;        /* result */
  331.  
  332.     if (qisfrac(q1) || qisneg(q1))
  333.         math_error("Non-positive integer required for qredcout");
  334.     rp = qfindredc(q2);
  335.     if (qiszero(q1))
  336.         return qlink(&_qzero_);
  337.     r = qalloc();
  338.     zredcdecode(rp, q1->num, &r->num);
  339.     if (zisunit(r->num)) {
  340.         qfree(r);
  341.         r = qlink(&_qone_);
  342.     }
  343.     return r;
  344. }
  345.  
  346.  
  347. /*
  348.  * Multiply two REDC format numbers together producing a REDC format result.
  349.  * This multiplication is done modulo the specified modulus.
  350.  */
  351. NUMBER *
  352. qredcmul(q1, q2, q3)
  353.     NUMBER *q1, *q2;    /* REDC numbers to be multiplied */
  354.     NUMBER *q3;        /* modulus */
  355. {
  356.     REDC *rp;        /* REDC information */
  357.     NUMBER *r;        /* result */
  358.  
  359.     if (qisfrac(q1) || qisneg(q1) || qisfrac(q2) || qisneg(q2))
  360.         math_error("Non-positive integers required for qredcmul");
  361.     rp = qfindredc(q3);
  362.     if (qiszero(q1) || qiszero(q2))
  363.         return qlink(&_qzero_);
  364.     r = qalloc();
  365.     zredcmul(rp, q1->num, q2->num, &r->num);
  366.     return r;
  367. }
  368.  
  369.  
  370. /*
  371.  * Square a REDC format number to produce a REDC format result.
  372.  * This squaring is done modulo the specified modulus.
  373.  */
  374. NUMBER *
  375. qredcsquare(q1, q2)
  376.     NUMBER *q1;        /* REDC number to be squared */
  377.     NUMBER *q2;        /* modulus */
  378. {
  379.     REDC *rp;        /* REDC information */
  380.     NUMBER *r;        /* result */
  381.  
  382.     if (qisfrac(q1) || qisneg(q1))
  383.         math_error("Non-positive integer required for qredcsquare");
  384.     rp = qfindredc(q2);
  385.     if (qiszero(q1))
  386.         return qlink(&_qzero_);
  387.     r = qalloc();
  388.     zredcsquare(rp, q1->num, &r->num);
  389.     return r;
  390. }
  391.  
  392.  
  393. /*
  394.  * Raise a REDC format number to the indicated power producing a REDC
  395.  * format result.  This is done modulo the specified modulus.  The
  396.  * power to be raised to is a normal number.
  397.  */
  398. NUMBER *
  399. qredcpower(q1, q2, q3)
  400.     NUMBER *q1;        /* REDC number to be raised */
  401.     NUMBER *q2;        /* power to be raised to */
  402.     NUMBER *q3;        /* modulus */
  403. {
  404.     REDC *rp;        /* REDC information */
  405.     NUMBER *r;        /* result */
  406.  
  407.     if (qisfrac(q1) || qisneg(q1) || qisfrac(q2) || qisneg(q2))
  408.         math_error("Non-positive integers required for qredcpower");
  409.     rp = qfindredc(q3);
  410.     if (qiszero(q1) || qisunit(q3))
  411.         return qlink(&_qzero_);
  412.     r = qalloc();
  413.     zredcpower(rp, q1->num, q2->num, &r->num);
  414.     return r;
  415. }
  416.  
  417.  
  418. /*
  419.  * Search for and return the REDC information for the specified number.
  420.  * The information is cached into a local table so that future calls
  421.  * for this information will be quick.  If the table fills up, then
  422.  * the oldest cached entry is reused.
  423.  */
  424. static REDC *
  425. qfindredc(q)
  426.     NUMBER *q;        /* modulus to find REDC information of */
  427. {
  428.     register REDC_CACHE *rcp;
  429.     REDC_CACHE *bestrcp;
  430.  
  431.     /*
  432.      * First try for an exact pointer match in the table.
  433.      */
  434.     for (rcp = redc_cache; rcp < &redc_cache[MAXREDC]; rcp++) {
  435.         if (q == rcp->num) {
  436.             rcp->age = ++redc_age;
  437.             return rcp->redc;
  438.         }
  439.     }
  440.  
  441.     /*
  442.      * Search the table again looking for a value which matches.
  443.      */
  444.     for (rcp = redc_cache; rcp < &redc_cache[MAXREDC]; rcp++) {
  445.         if (rcp->age && (qcmp(q, rcp->num) == 0)) {
  446.             rcp->age = ++redc_age;
  447.             return rcp->redc;
  448.         }
  449.     }
  450.  
  451.     /*
  452.      * Must invalidate an existing entry in the table.
  453.      * Find the oldest (or first unused) entry.
  454.      * But first make sure the modulus will be reasonable.
  455.      */
  456.     if (qisfrac(q) || qiseven(q) || qisneg(q))
  457.         math_error("REDC modulus must be positive odd integer");
  458.  
  459.     bestrcp = NULL;
  460.     for (rcp = redc_cache; rcp < &redc_cache[MAXREDC]; rcp++) {
  461.         if ((bestrcp == NULL) || (rcp->age < bestrcp->age))
  462.             bestrcp = rcp;
  463.     }
  464.  
  465.     /*
  466.      * Found the best entry.
  467.      * Free the old information for the entry if necessary,
  468.      * then initialize it.
  469.      */
  470.     rcp = bestrcp;
  471.     if (rcp->age) {
  472.         rcp->age = 0;
  473.         qfree(rcp->num);
  474.         zredcfree(rcp->redc);
  475.     }
  476.  
  477.     rcp->redc = zredcalloc(q->num);
  478.     rcp->num = qlink(q);
  479.     rcp->age = ++redc_age;
  480.     return rcp->redc;
  481. }
  482.  
  483. /* END CODE */
  484.